Measuring the fluctuation-dissipation ratio in glassy systems with no perturbing field 
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A method is presented for measuring the integrated response in Ising spin system without applying 
any perturbing field. Large-scale simulations are performed in order to show how the method 
works. Very precise measurements of the fluctuation-dissipation ratio are presented for 3 different 
Ising models: the 2-dimensional ferromagnetic model, the mean-field diluted 3-spin model, and the 
3-dimensional Edwards- Anderson model. 



Disordered and frustrated models are a fascinating but 
still poorly understood subject in contemporary statisti- 
cal mechanics. The interest in these systems also comes 
from their many interdisciplinary applications: from the 
physics of glass-former liquids to that of polymers and 
biomolecules, from the description of error correcting 
codes to the study of the computational complexity and 
phase transitions in theoretical computer science. 

Here we will use the term glassy system for a generic 
model showing very slow relaxation to equilibrium 
Because of the huge equilibration time, a glassy system 
may be in the out of equilibrium regime for all the ex- 
perimental times. Then a complete understanding of this 
regime is what one needs in order to correctly describe a 
real slow-evolving material. Moreover, numerical studies 
of the off-equilibrium regime do not suffer from finite-size 
effect since very large sizes can be used. They present fi- 
nite time corrections which can be usually kept under 
control, thus allowing for better numerical estimations. 

Among the numerical methods that can be used in 
the out of equilibrium regime, the study of the so-called 
fluctuation- dissipation ratio (FDR) 0] has been shown to 
be a very successful one 0, 0] . This method is based on 
the comparison of how spontaneous and induced fluctu- 
ations relax. Actually one measures an autocorrelation 
function C(t, s) 01 an d the associated response function 
R(t, s) and defines the FDR X(t, s) through the formula 

TR(t,s) = X{t,s)d 3 C(t,s) , (I) 

where T is the temperature. At equilibrium the fluctua- 
tion dissipation theorem (FDT) holds, implying X = 1. 
In the large times limit — s, t — > oo with C(t, s) — > q 
- the FDR X(t,s) converges to the limiting function 
X(q). The physical meaning of the function X(q) has 
been explained in Refs. jjij], where it has been shown that 
under some hypothesis (stochastic stability) the equation 

X(q)= X (q)= f ' P(q')dq' (2) 
Jo 

holds. In Eq.(J2J P(q) represents the overlap pdf in the 
threshold states, that is the states reached by the out of 
equilibrium dynamics on very large times, which could be 



different from the thermodynamical state p| . It has been 
conjectures that the effective temperature T e g = T IX 
plays a central role in off-equilibrium glassy systems [6(. 

In numerical simulations the punctual response func- 
tion i?(i, s) is very noisy, while a much better signal can 
be obtained for the integrated response function 

X (t,t w )=T^ R{t,s)ds . (3) 

With respect to the usual definition, the temperature T 
has been added in the above equation in order to simplify 
the notation in the following formulae and to have a well 
defined x(t,t w ) in the T — > limit. In the large time 
limit, Eq.lO can be rewritten as 

X (C)= f X(q)dq . (4) 

Jc 

So the FDR can be simply written as X(C) — —dcx(C)- 
The aim of this Letter is to propose and to show the 
efficacy of a very precise method for measuring the inte- 
grated response x(C) and the FDR X{C) in spin models. 

Up to now the best protocol for measuring x(C) in 
spin systems has been the following 0,0: 

• initialize the system in a random configuration; 

• quench the system at a temperature T < T c and 
evolve it for t w Monte Carlo sweeps (MCS); 

• switch on a random magnetic field of very small 
intensity h and continue evolving the system while 
measuring C h (t,t w ) and Xh(t,t w ). 

The parametric plot of Xh (t, t to) versus Ch,(t>tw*) con- 
verges to the function x(C) in the limit t w — > oo and 
h — > 0. Even when extrapolations can be safely done, 
they always require a large numerical effort: for example, 
in order to correctly take the h — ► limit, the whole sim- 
ulation must be repeated for many h values in the linear 
response regime. Moreover in frustrated systems like spin 
glasses the response may have strong non-linearities even 
for very small probing fields and it is usually very hard 
to predict a priori which is the linear response regime. 
Furthermore in out of equilibrium simulations the size 
of the linear response regime may change with the age 
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of the system: A fair conjecture is that it decreases for 
larger times. If this would be true, extrapolations to the 
interesting limit would become still more difficult. 

For all these reasons we consider of primary impor- 
tance the development of a method which allows one to 
calculate the linear response in a spin system without ap- 
plying any probing field. After having taken analytically 
the h — -> limit, one is left only with the t w — » oo limit. 
This limit will be somehow unavoidable as long as the 
only way for aging a glassy system will be to simulate it 
for a long time |l3| . 

Inspired by a recent work by Chatelain we write 
down an analytical expressions giving the integrated re- 
sponse x(t, tw) in a simulation with no probing field [l4j . 

Let us specialize on systems with N Ising spins and 
Hamiltonian H.q (generalization to Potts variables is 
straightforward 8J). The Hamiltonian 7io may contain 
some quenched disorder, but we do not need to specify 
it, since our calculations hold for a generic TIq, either 
disorder or not disordered. In the former case the final 
result can be eventually averaged over the quenched dis- 
order distribution, but following formulae are valid for 
any given disorder realization. 

When the probing field is switch on the Hamiltonian 
becomes Ti = Tig — ^i cr ij where hi are i.i.d. random 
variables with hi — and hihj — h 2 Sij. For simplicity 
we define hi — hei with £7=0 and el e] = Sij . 

The FDR for the observable Ait) — ^ i Ei<Ji(t) is given 
in terms of the correlation and response functions 



NC(t, s) 
NR(t, s) 



(A(t)A0O)=X>(tM5)), 



(5) 



d(A(t)) 
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d((Ti(t)) dhj 



dhj(s) dh 
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dhjis) ^ dhi(s) 



(6) 



where ( • ) represents the average over thermal histories. 

We use a discrete-time dynamics as the one taking 
place in a Monte Carlo simulation. The time t counts 
the number of single spin updates and not the number of 



Monte Carlo sweeps (which is then t/N). The function 
I(t) gives the index of the spin to be updated at time 
t, and so it depends on the updating rule (e.g. random 
or sequential). At the i-th time step the spin Oi with 
i = I(t) is updated according to heat-bath probabilities 



prob {<Ji 



exp^fef + hj)] 
2cosh[/?(hf +hi)] 



(7) 



where j3 is the inverse temperature and the Weiss field 
hj^ takes into account the effect of Hamiltonian Ho on 
the spin to be updated. For example, in the case of 2- 
spin interacting Hamiltonians the Weiss field is given by 



hf 



Under this dynamics the expectation value of the j-th 
spin a time t is given by 



(<Tj(t)) = Ti 3(tl) 



<rj®tlW m (ff(lf)\S(!?-l) 



(8) 



where a is a short-hand notation for the N spins config- 
uration, the trace is over all the trajectories S(t') with 
1 < if < t, and the transition probability is given by 



Wi{ff\r) 



eMP^(hT + hi)} 
2cosh[/3(/if + hi)} 



(9) 



Note that hj*(cr) = hj*(r) since it does not depend on 
the spin in i. The transition probability Wi only depends 
on the perturbing field on site i, such that 



dWi{a\f) 



dhj 



kj W t {a\f)f3 (at 



(10) 



h=0 



where we have defined = tanh(/3/iJ v ). 

Now we suppose that an infinitesimal probing field hk 
on site k is switch on at time t w : hk(t) = h9(t — t w ). 
This means that the transition probability Wk (and only 
this one) will depend on the perturbing field for all times 
larger than t w . Differentiation of Eq. © with respect to 
this field yields the integrated response 



Xjk(t,t w ) = T 



d(a 3 (t)) 



h=0 



(11) 



t' = l s=t w + l 

t 

<Tj(t) Aa k (t,t w )) with Aa k (t,t w )= s i(s),k (°k(s) ~ ^{s)) . 

r 



The correlation in Ea. Hllfl is what one has to measure in to get the integrated linear response. 

a numerical simulation with no perturbing field in order Few comments are in order. The time-integrated quan- 
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FIG. 1: FD plot for the 2D Ising ferromagnet at T = 2. The 
horizontal line is the analytic long time limit. 



FIG. 2: FD plot for the long-range 3-spin model with fixed 
connectivity 4 at T = 0.5. The line 0.53745 - 0.50256 C is the 
best linear fit to t w — 10 3 data with C < 0.9. 



tity AcTfe only gets contributions when the spin <Jk is up- 
dated, so most of the times it is unchanged. The con- 
tributions summed up in Aofc are the differences among 
the actual value of the spin ak and the expected one 
aW . So Act/c is a random variable with zero mean, 
(<7fc) = a™ => (A<7fc) = 0, and variance (Act 2 ) oc (t — t w ). 

In the T — limit Ea.pijl has a nice and simple 
physical interpretation. Since for T = we have that 
a k ~ = a k^h w .0' then Acfe takes a contribution 
only when the spin <7k has a zero Weiss field on it, i.e. 
it is free to respond to an infinitesimal field. If the 
Weiss field is different from zero the spin is completely 
frozen and it can not respond to an infinitesimal per- 
turbing field. So the integrated response in Ea. l|ll|l can 
be rewritten as a simple sum of correlation functions 
Xjk{t,t w ) = (t)<Tfc (s)} , where the primed sum is 

over all the times larger then t w when au is updated un- 
der a zero Weiss field, i.e. being a free spin. 

We now present numerical results for 3 spin models 
which are believed to belong to 3 different classes: ferro- 
magnetic model in 2D (coarsening system), diluted long- 
range 3-spin model with fixed connectivity 4 (discontinu- 
ous spin glass) and Edwards- Anderson model in 3D (con- 
tinuous spin glass). For each model we have checked 
that the Xhi^h) curve measured with the perturbing field 
converges for h — ► to the one measured with the new 
method. Hereafter times will be counted in MCS. 

The first model is the ferromagnetic Ising model on 
the 2-dimensional square lattice. We have simulated at 
T = 2 ~ 0.88 T c systems of sizes 1000 2 and 7000 2 in order 
to check the absence of any finite-size effect (the data we 
show are from the 1000 2 samples). For each waiting time, 
t w = 10 2 ,10 3 ,10 4 , averages have been taken over 100 
different thermal histories, and the corresponding x(C) 
curves are shown in Fig. ^ The horizontal line is the 
analytical prediction for the large times limit, \ = 1 — 
m cq = 0.17. Numerical curves are clearly compatible 
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FIG. 3: FD plot for the 3D Edwards-Anderson model. 



with the analytical asymptote in the large times limit. 

The second model we studied is the 3-spin model de- 
fined on a random hypergraph with fixed connectivity 4. 
This model has been solved analytically with a one-step 
replica symmetry breaking Ansatz in Ref. [ll| . The dy- 
namical critical temperature is Td — 0.755 ± 0.01. We 
have run simulations for a size N = 999999 at tem- 
perature T — 0.5 ~ 0.66 and the resulting x(C) 
curve is shown in Fig. [21 The number of samples used 
is 10 for t w = 10, 10 2 , 50 for t w = 10 3 and 20 for 
t w = 10 4 . The straight line in Fig. [21 is a linear fit to 
t w = 10 3 data in the region C(t,t w ) < 0.9, which per- 
fectly interpolates the data (x 2 /d.o.f. = 0.82). It gives 
a Parisi breaking parameter on threshold states equal to 
mthiT = 0.5) = 0.5 ± 0.02. The error is an estimate 
of systematic effects, mainly given by the slight increase 
of m with t w . Comparison of this value for m t h with 
corresponding static predictions will be done in Ref. |lll | . 

The third model we studied is the 3-dimensional 
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FIG. 4: Same as Fig. El with rescaled variables. Inset: FDR 
for T = 0.75 obtained from the derivative of the rescaled data. 



Edwards- Anderson model with J = ±1 couplings, which 
undergoes a phase transition to a spin glass phase at 
T c = 1.14 ± 0.01 [uj. We have simulated samples of 
size L = 20 at temperatures T = 0.75 ~ 0.66 T c and 
T = 0.5 ~ 0.44 T c , for three different waiting times 
t w = 10 2 , 10 3 , 10 4 . The results are shown in Fig. El 

For a given temperature the x(C) curves look very sim- 
ilar in shape, the main difference being the ^-dependent 
Edwards- Anderson order parameter q EA {t w ), here de- 
fined as the point where the x(C) curve leaves the 
FDT line 1 — C. In order to exploit all the data we 
tried to collapse the curves before fitting. The collapse 
can be achieved either by shifting the curves such that 
the q EA (t w ) coincide, either by the following rescaling: 

C Ie s(t,tw) = ^C(t,t w )/q EA (t w ),Xrcs(t,t w ) = 1 — A (1 — 

X(t,t w ))/q EA (t w ), with an arbitrary A. Both scaling are 
statistically acceptable. In Fig. 0] we show the second one 
which is slightly better, with A = q EA (10 4 ). 

If the measured data are already in the asymptotic 
regime, i.e. the scaling is valid for larger times, and 
since lim^ — >ooQ EA (t w ) = <7ea > 0, we can conclude that 
the FDR is non trivial in the 3-dimensional Edwards- 
Anderson model, with an X(C) like the one depicted in 
the inset of Fig. H for T = 0.75. 

The Edwards-Anderson model is the one which took 
the great part of the simulation time. Indeed, in order 
to have reasonable error bars, we ran at each of the two 
temperatures 10 4 samples for t w < 10 3 and almost 3 • 10 4 
samples for t w = 10 4 . Solely the t w — 10 4 runs took the 
equivalent of more than 1 year of CPU-time on a latest 
generation 2.0 GHz computer. This is a consequence of 
the fact that errors on the linear response x(t,i UJ ) in- 
creases like y/t — t w and so the number of samples for 
keeping the error on x(C) constant increases more or 
less linearly with the waiting time t w . We believe this is 
the main drawback of the present new method for mea- 
suring the linear response: Although it is very successful 



for small times, it becomes very noisy at larger times and 
so it requires a huge statistics. 

From this observation one could conclude that the 
usual old method of measuring the response with a small 
perturbing field would eventually remain the only valid 
one, but this is not true. Very probably the linear regime 
in the perturbing field h decreases with the age of the 
system. In order to remain in the linear response regime 
one should decreases the intensity of the perturbing field 
during the simulation, thus increasing the error on the \ 
for late times. 

Let us conclude with two remarks. First remark: Hav- 
ing understood that the integrated response can be writ- 
ten as a correlation function, it should be clear that all 
the functions C(t,t w ) and x{tjtw) can be calculated in 
the same simulation for any value of t and t w . Moreover, 
being correlation functions self-averaging quantities, it 
should be possible in principle to calculate them in a 
single simulation of a sufficiently large sample. Second 
remark: The method presented here can be also used 
for any other Monte Carlo simulation (e.g. glass-former 
particle systems). The only condition for using present 
analytical expressions is the discreteness of time. 
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